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Abstract. We study the dynamics of a spherical rigid body that rocks and 
rolls on a plane under the effect of gravity. The distribution of mass is non- 
uniform and the centre of mass does not coincide with the geometric centre. 
The symmetric case, with moments of inertia I\ = I2 < I3, is integrable and 
the motion is completely regular. Three known conservation laws are the total 
energy E, Jellett's quantity Qj and Routh's quantity Qr. When the inertial 
symmetry I\ = I2 is broken, even slightly, the character of the solutions is 
profoundly changed and new types of motion become possible. We derive the 
equations governing the general motion and present analytical and numerical 
evidence of the recession, or reversal of precession, that has been observed in 
physical experiments. We present an analysis of recession in terms of critical lines 
dividing the (Qr,Qj) plane into four dynamically disjoint zones. We prove that 
recession implies the lack of conservation of Jellett's and Routh's quantities, by 
identifying individual reversals as crossings of the orbit (Q n(t) , Q j (t)) through 
the critical lines. Consequently, a method is found to produce a large number of 
initial conditions so that the system will exhibit recession. 



1. Introduction 

We investigate the dynamics of a spherical rigid body rolling on a plane. The 
distribution of mass is non-uniform, so that the centre of mass does not coincide 
with the geometric centre. However, the line joining the mass centre and geometric 
centre is assumed to be a principal axis. We denote the principal moments of inertia 
by Ii, I 2 and ^3, and assume that I\ < I 2 < I3. The symmetric case, when I\ = I2, 
was first studied by Routh [15], and in this case the body is called Routh's Sphere. 
There are three constants of motion and the system is integrable. In the asymmetric 
case, 1\ 7^ I2, the system is no longer integrable. We find that even a small degree of 
asymmetry has a dramatic effect on the motion of the body. 

The equations of the symmetric loaded sphere are identical to those governing the 
motion of the tippe-top, which has been studied extensively (see [7] for a comprehensive 
reference list). However, in the case of the tippe-top, the angular momentum about 
the principal axis with maximum moment of inertia is large, and sliding friction plays 
a key role. In the case under consideration here, we are interested in solutions where 
the angular velocity remains moderate and there is pure rolling contact. There are 
two characteristic modes of behaviour: pure rocking motion in a vertical plane, and 
pure circular rolling motion. The general motion has aspects of both these special 
cases, which leads us to name the body the rock'n'roller. 
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Figure 1. The physical rock'n'roller, constructed by slicing off a polar cap from 
a standard bowling ball. The polar angle is © ?s 53° . 



This investigation arose from the observation of the oscillations of a glass candle 
holder, spherical in form with an opening at the top. For a more systematic study, we 
constructed a larger and more massive body by removing a polar cap from a bowling 
ball to produce a truncated sphere (figure 1). As long as the tilting angle is such that 
the geometric centre is vertically above the contact point, the dynamics are equivalent 
to those of a loaded sphere. It was found that when the ball was tilted over to an angle 
of about 130°, it rocked back and forth but also precessed through an azimuthal angle 
that alternately increased and decreased. This unexpected and surprising recession, 
or reversal of precession, demanded an explanation in terms of dynamics. 

We will show that for a symmetric loaded sphere reversal of the precession is 
impossible. This raises the question: what factor is missing from our dynamical 
model? We rule out sliding friction, since the motion is gentle with no evidence of 
slipping. Random perturbations, due to the imperfect shape of the ball or irregularities 
of the underlying surface, were not considered as a likely cause of the behaviour, as 
experiments indicated that the recession was quite a robust feature of the motion. 

Although bowling balls are manufactured to high tolerence, and deviations from 
perfect sphericity must be very small, slight anomalies in the mass distribution are 
unavoidable. Moreover, the recesses in the physical body remaining from the finger 
holes introduce some asymmetry (figure 1). We were thus led to study the dynamics 
when the inertial symmetry I\ — I2 is broken. We find that even a minute deviation 
from symmetry changes the behaviour of the numerical solution profoundly. Of the 
three quantities conserved in the symmetric case (total energy E, Jellett's quantity 
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Figure 2. Hierarchy of loaded spheres. The vector Of) is from the mass centre 
O to the geometric centre C, and k is the unit vector along the /3-axis. See text 
for details. 



Qj and Routh's quantity Qr), only the energy remains invariant when I\ ^ 1%. 
We derive the equations governing the general motion and present analytical and 
numerical evidence of recession. We base our analysis on the existence of critical 
lines dividing the (Qr, Qj) plane into four dynamically disjoint zones. We prove 
that recession implies the lack of conservation of Jellett's and Routh's quantities, 
by identifying individual reversals as crossings of the orbit (Q^(t), Qj(t)) through the 
critical lines. This leads to a method of defining initial conditions for which the system 
will exhibit recession. 

The rock'n'roller is one of a hierarchy of loaded spheres. For the most general 
case, the vector Od from the mass centre O to the geometric centre C does not lie 
on a principal axis, and all moments of inertia are distinct. This is called Chaplygin's 
Top [3]. For the rock'n'roller, the geometric centre lies on a principal axis and 06 
is parallel to k, the unit vector along the /3-axis. Routh's Sphere is the special case 
of this with 1\ = I2 and Chaplygin's Sphere the special case where the mass centre 
and geometric centre coincide. The hierarchy is illustrated in figure 2. For recent 
discussions, see [9, 5, 4, 10, 17, 2, 16] 

See an animation of the rock'n'roller that exhibits precession and recession in 
http://mathsci.ucd.ie/~plynch/RnR/RnR_movie.gif. We produced this movie 
from a Mathematica simulation code of the equations in the asymmetric case I\ 7^ 1% , 
corresponding to the initial conditions described in Figure 8. 
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Figure 3. The intermediate coordinate frame used to study Routh's Sphere. 



2. Symmetric Body (I ± = J 2 ): The Dynamical Equations 

We consider a body, spherical in shape with unit mass and unit radius, whose mass 
distribution is non-uniform but symmetric about some line through the centre. We 
assume that the centre of mass is off-set a distance a from the geometric centre and 
that the moments of inertia perpendicular to and along the symmetry axis are I\ and 
Is, with h{= h) < h- All the parameters are determined once the angle of the polar 
cap that is removed is known (see Appendix §A.l). In an inertial frame of reference, 
the equations governing the dynamics of the body are 

df v ; 

where v is the velocity of the centre of mass in the absolute frame and F the total 
force acting on the body; and 

where L is the intrinsic angular momentum and G the total moment about the centre 
of mass. 

The derivation in this section is similar to that in [19]. We consider a rotating 
frame of reference, with unit triad (i',j',k') whose origin moves with the centre of 
mass of the body. The vector k' is aligned with the axis of symmetry of the body 
and j' is in the same vertical plane as k' (see figure 3). Then i' is horizontal and 
perpendicular to the plane of the figure, pointing inward. We use primes for this 
intermediate frame to distinguish it from the body frame that will be introduced in 
§4 below. 
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The angular velocity of the body, expressed in the intermediate frame, is 

u) = + ui' 2 j' + u' 3 k' 

Although this frame is not fixed in the body, it forms a set of principal axes in the 
symmetric case, I\ = h, and the angular momentum is given by 

L = huj'^i' + hw 2 j' + huJ 3 k' . 

We denote the angular velocity of the frame itself by fi and note that 

n = 9i' + 4>K = 9i' + s<j>j' + c4>k' = (Qi,Q 2 ,0 3 ) ( 3 ) 

where s = sin 9 and c = cos# and K is a unit vertical vector. The Euler angles 
(9, <f>, tjj) are related to the components of angular velocity by 

u[ = 9 , uj' 2 = s<j) , u}' 3 ~c(j) + ip. (4) 

Definitions are standard, and may be found in [12, 18, 20]. For a list of the principal 
symbols used in this study, see Table 1. 

2.1. Equations in the intermediate frame 

In the moving frame, the equations (1) and (2) become 

^- + nxv = F (5) 

at 

and 

^ + OxL = G (6) 
Expanding these in components in the i'j'fe'-frame, we get 

v 2 + tt 3 v[ - Qxv's = F 2 (7) 

v' 3 + niv' 2 - n 2 v[ = f 3 

for momentum. The angular momentum equations become 
I\u)[ + h^l2^' 3 - = G\ 

I 2 lu' 2 + Ji0 3 wi - h^li^ = G 2 (8) 
hu' 3 = G 3 

Eqns. (7) and (8) are identical to (12.412) in [18] (with h = h). 

The forces acting on the body are gravity W — (0, — gs, —gc) and the force of 
reaction R = (i?i , R 2 , R 3 ) : 

F = W + R 

Defining / = c — a, the vector from the point of contact P to the centre of mass 
O is r = (0, s,f) (see figure 3). Then the total moment about O is given by 
G = r X R. The constraint of no slipping at the contact point requires that the 
body is instantaneously rotating about this point. Thus, 

v = w X r = {fJ 2 - su' 3 , -fu[, su}[) . (9) 

The reactive forces may be eliminated by combining the angular momentum equation 
(6) with the vector product of r and the momentum equation (5). The velocity v 
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Table 1. Principal symbols used in this study 



G 



Symbol Meaning 

F Total forcing in Newton's equation 

G Total moment in Newton's equation 

h i h, h Principal moments of inertia of body 

K Unit vertical vector 

L Angular momentum of body about centre of mass 

C Lagrangian function 

Qj Jcllctt's quantity, constant in symmetric case 

Qr Routh's quantity, constant in symmetric case 

R Force of reaction at contact point 

T Total kinetic energy 

V Potential energy 

V Velocity of centre of mass in the space frame 
W Force due to gravity (weight) 

a Distance from geometric centre to centre of mass 

c Cosine of tilting angle, c = cos 9 

Cosine of azimuthal angle, = cos <p 

d Cosine of polar angle, d = cos 

/ Projection of vertical radius on fc-axis, / = cos 9 — a 

g Acceleration of gravity 

h Height of centre of mass, h = 1 — a cos 9 

i,j,k Principal unit orthogonal triad in body coordinates 

k' Principal unit orthogonal triad in body coordinates 

r Moment vector, from contact point to mass centre 

s Sine of tilting angle, s = sin 9 

Sine of azimuthal angle, s^ = sin <fi 

t Time 

v Velocity of centre of mass in the body frame 

Vi, V2, V3 Components of v in body coordinates 

Co- latitude of polar cap removed to construct the rock'n'roller. 

T Rotation matrix 

$ Azimuthal angle spanned by solution, $ = max — 4> m i n 

CI Angular velocity of intermediate frame 

e Asymmetry parameter, e — (I2 — 

9, <f>, ip Euler angles (tilting, azimuth and spinning angles) 

fik Lagrange multipliers 

p Measure quantity, p = [I 3 + s 2 + {h/ h)P]~ 1/2 

a Sine of spining angle, a = sin-0 

t Period of rocking motion 

X Cosine of spinning angle, \ = cos ip 

u> Angular velocity of body 

u[, ll> 2 , ^'3 Components of u> in intermediate coordinates 

u)\, ll>2, 0J3 Components of u; in body coordinates 
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may be expressed in terms of the rotation u> by means of the constraint (9). We then 
obtain three equations for Wi, uo 2 an d ^3: 

'/i + .s 2 + / 2 







-fa 
h + s 2 











- 








V^3 



(10) 





h + P 
-fs 13 

where Pi, P 2 and P3 depend on the angles and angular velocities. Full details of the 
derivation are presented in the Appendix (§A.2). The rates of change of the angular 
variables follow from (4): 

u 2 /s, ip = uj' 3 - {c/s)uj 2 . (11) 



= w. 



We now have six equations (10) and (11) for the six variables {9, (p, ip, co[, u> 2 , W3}. 



2.2. Special solutions 

Pure Rocking. For pure rocking motion, with no change of azimuthal angle and no 
rotation about the axis of symmetry, we have <j> = ip = and so u> 2 = lo' 3 
the system reduces to a single equation for the tilting angle 9: 



0. Then 



(g + o 2 )a 



sin 6 = 



h + f 2 + s 2 

For small amplitude 9 <^ 1, and assuming a <C 1, this becomes 
ga 



= 0. 



./1 + 1_ 

the equation for simple harmonic oscillations. 



(12) 



(13) 



Pure Rolling. For the case of pure circular rolling motion we have 

9 = 0, (f) = constant , ip = constant 

so that Oi = u>[ = and io' 2 and cu' 3 are constants. It follows immediately that 
P 2 = P 3 = (sec (A. 2)). The requirement that 9 = 6 , constant, implies Pi = 0, 
which yields a relationship between u 2 and uj' 3 : 

, _ (h cot 9 + mh f esc 9 )uj 2 2 - ga sin 9 

where /o = cos6*o — a and ho = 1 — acos#o are constants. If we start with an d 
u 3 related by (14) and 9 slightly perturbed from 9q, motion with nutation about 9o 
results. 



2. 3. Constants of motion and general solution 

We consider the case of a perfectly rough contact, with rolling motion. Given that 
there are two symmetries in the problem, invariance under addition of arbitrary 
constants to either <j> or ip, we might expect two invariants in addition to the total 
energy. For general initial conditions, there are three constants of integration. They 
arc the total energy, Jellett's constant and Routh's constant (see [7] for a complete 
derivation of these constants). 

The kinetic energy is the sum of translational and rotational components: 

T = ±[< + v 12 + v' 2 } + \[hJ 2 + IW 2 + I 3 cu' 2 } 
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and the potential energy is 

V = ga{l-cos6). 
Then, since there is no dissipation, the total energy 

E = T + V. (15) 

is conserved. Jellett's constant is the scalar product of the angular momentum and 
the vector joining the point of contact to the centre of mass: 

Qj = L-r = I 1 8u' 2 + I 3 fu' 3 , (16) 

and Routh's constant, more difficult to interpret physically, is 

Qn = f (17) 
where, following [2], we define the measure 

p{6) = l/V/3 + a 2 + (/ 3 //i)/ 2 . (18) 

Notice that our definition of Routh's constant differs from the usual quadratic function 
of u> 3 , in [15], [7] and elsewhere, 

An interesting historical discussion of these constants may be found in [7]. Note 
that the constancy of Qr implies conservation of the sign of u>' 3 : since the measure p 
is positive definite, u)' 3 cannot pass through zero. For the tippe-top, this precludes the 
tipping phenomenon for the case of rolling motion. 

From the equations (11) determining the rates of change of the angles, we can 
solve explicitly for <fi and ip in terms of 9 and Jellett's and Routh's constants: 



= V^6,Qj,Q R ) = -^ 

1\S 



Qj - PfhQR 



^ = V^0,Qj,Q r ) = -j^ 

1\S 



cQj-p(cfh + hs 2 )Q R 



(19) 



(20) 



Since Qj and Qr arc constants, the rates of change <fi and ip are determined as single- 
valued functions of the angle 9. We will show in the next section that recession, 
or reversal of precession, implies in particular that <j> and tp at a given angle 9 
systematically change their sign as time evolves. Therefore, in the symmetric case 
it is impossible to have recession for Routh's Sphere. 

We can use the constants of motion to reduce the system to a single equation 
for the tilting angle 9. We use Routh's constant (17) to obtain lo' 3 (9). Then Jellett's 
constant (16) gives uj' 2 (9). Finally, the energy (15) gives an expression for ^[(9), 
yielding an equation of the form 

9 2 = F{9), (21) 

which may be integrated to obtain 9(t). As a result, the system can be explicitly 
integrated. However, we will not derive explicit analytical expressions for F{9) and 
9(t). The reader is referred to [7] for a more explicit treatment; see also [4], [16]. We 
see that the evolution of 9(t) obtained from (21) gives the rocking component of the 
motion, while the evolution of 4>{t) and obtained from (19) and (20), give the 

rolling and spinning components of the motion. 
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2.4- Precession of the Rocking Motion 

The generic motion of the symmetric body is quasi-periodic. On the one hand there 
is a period t of the rocking motion, determined by the equation 

6x d6> 

where < 9n <tt, and 9^ and 9x are the turning points where F(9n) = and 

F(9x) = 0. On the other hand, the rolling motion during this period can be computed 
by integrating the rates of change of the angles cf) and ip from (19) and (20): 

A0 = f V+W), Qj, Q R ) At = 4 f X Me ^J t a? R) dd ' ( 22 ) 
Jo Je N \JF(v) 

with an analogous formula for the angle A-0. Gcncrically, A<p is not commensurate 
with 2tt; this implies the quasi-periodicity of the precessing motion. As a consequence 
of quasi-periodicity, the projection of the trajectory onto the #-<j()-plane densely covers 
a two-dimensional region. 

In order to quantify the precession, we distinguish two angles, the full azimuthal 
angle and the visible angle </>(mod27r), which is the angle that is seen by an 
observer in the space frame. (We will occasionally use the visible half-angle <p (mod tt), 
which gives more illustrative plots in the case of the (asymmetric) rock'n'roller). 
Correspondingly, there will be two types of precession angle: A</>, the full precession 
angle defined by (22), and A(f> (mod 2tt) , the visible precession angle. 

2. 5. Qualitative analysis of the precession. Criticality 

We now estimate the precession angles A0 and A0(mod27r) from (22). Heuristically, 
the main contribution comes from the regions near the turning points, where F(9) = 0. 
The relative contributions at 9^ and 9 X will be determined by the magnitude and sign 
of V(f,(0 , Q j , Q r) at the turning points. It is therefore useful to study separately the 
behaviour of V^(9 7 Qj, Qr) for 9 in each of the asymptotic regions 9 w and 9 w 7r: 

• Asymptotic region 9 w 0. A Laurent expansion of (19) gives 

W, Qj, Qr)\ 6 ~o = @J Q"fi) + °« > (23) 
where we define the 'critical Jellett quantity at 9 = 0' as: 

Q<y^ = p (l-a)I 3 Q R , (24) 

with p = l/y/I 3 + (h/h)(l-a) 2 . 

• Asymptotic region 9 w ir. A Laurent expansion of (19) gives 

V,(0,Qj,QR)\ g ^ = h ^_ ey (Qj-Qt)+0{1), (25) 

where we define the 'critical Jellett quantity at 9 = 7r' as: 

Q c j^ = -p4l + a)I 3 QR. (26) 

with P7r = \/y/h + {h/h){\ + a) 2 . 

• Monotonicity property. The factor [Qj — pf I 3 Qr] appearing on the right-hand 
side of (19) is a monotonic function of the angle 9 £ [0, n]. The proof of this is 
straightforward. 
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From the above asymptotic expansions, we conclude that, in the space of initial 
conditions parameterised by (Qr,Qj), there are four regions of interest, and the 
behaviour of 4> is qualitatively different in each region: 

Region I: Q R > 0, Qj£ < Qj < Q C J$. The function V<f,(0, Qj, Qr) goes from -co 

at 9 = to oo at 9 = n. From the monotonicity property it follows that this 

function has a single zero. 
Region II: QJ* < Qj, QJfi < Qj. The function V c f > (9,Qj,QR) goes from oo at 

9 = to oo at 9 = 7r. From the monotonicity property it is possible to show that 

this function is positive definite and has a single minimum. 
Region III: Q R < 0, QJ* < Qj < QJ*. The function V^.Qj^r) goes from oo 

at 9 — to — oo at 9 = n. From the monotonicity property it follows that this 

function has a single zero. 
Region IV: Qj < QJ*, Qj < QJ*. The function V^O, Qj, Qr) goes from -oo at 

9 = to — oo at 6 = it. From the monotonicity property it is possible to show 

that this function is negative definite and has a single maximum. 

Similar results can be obtained for the velocities tp, but these are omitted here. 

In figure 4 we show the four regions, separated by the two critical lines Qj = QJ* 
(solid red line from top left to bottom right) and Qj — QJ* (dashed blue line, from 
bottom left to top right). Typical plots of the function V < p(9,Qj,Qr) versus 9 are 
inserted in each region. The asymptotic behaviours are evident. 

The critical Jellett quantity QJ* plays a key role in determining the visible 
precession angle A(/>(mod27r) in the interesting case Ox ~ tt- The main contribution 
to the precession angle comes from the turning point 9 = 9x and from (25) we see 
that the sign of this contribution depends on the sign of Qj — QJ*- For example, 
if an initial condition with Qj — QJ 1 * <; (Regions I or II) has a precession angle 
A0 = ao (mod27r), then a slightly different initial condition with Qj — QJ 1 * ~ 
(Regions IV or III) will have a precession angle Acj) — — a (mod 2n); the corresponding 
motion will appear to be reversed. 

The critical Jellett quantity QJq determines the full precession angle A(f> when 
9 jy w 0. The main contribution to A<p comes from the turning point 9n, and is 
given by rapid changes of <j> in jumps of approximately ±7r, the sign of these jumps 
depending on the sign of Qj — QJq- In this way, an initial condition in Region I or 
IV will give rise to a full-angle precession A(f> < whereas, for initial conditions in 
Region II or III, A</> > 0. It is worth mentioning that this critical quantity is related 
to the energy of the system since it appears in the Laurent expansion of the function 
F(9) near 9 = 0. Sec [7], where this critical quantity was identified in terms of the 
centrifugal barrier. 

2. 6. Quantitative estimate of applicability of criticality criteria 

Let us consider the asymptotic region 9 ~ n. For the above asymptotic analysis to 
be of practical importance, the maximum rocking angle 6x must be close to 7r. Only 
then will the asymptotic Laurent expansion (25) determine, to a good approximation, 
the value of <j> at = Ox- In particular, we will observe a dramatic difference in <j) at 
9 = 9x and in the precession angle when considering two nearby points, one in Region 
I and one in Region IV. 

To quantify how close should Ox be to ir, necessary conditions are: (i) In Regions 
I and III, 6> z , the zero of <j>{9), must be less than 9x] (ii) In Regions II and IV, 9 C , 
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Figure 4. (Colour online). The four critical regions defined by the critical lines 
Qj = Q C J^ (red, top left to bottom right) and Qj = Q c Jq (blue, bottom left to 

top right). In each region, the graph of as function of 6 is shown, for selected 
values Qr, Qj and parameters a = 0.05, ^3 = 2/5 and Ii = (1 — 5 a/2)Lj,. 



the extremum of (f>{9), must be less than Ox- In each case, there is a relation between 
&x, Qr, Qj and the parameters a, I\, I3. 

For Regions I and III, this condition has a simple analytical formulation: 



, . , . , aPh + Vhjh - Wgh - h)(h + 1) ~ a 2 Ph + hh(h + 1) 
-1 < cos Ox < cos6» z = — , 

P\h -h) + hh 

where (3 = I3 — (Qj/Qr) 2 - Realistic values of parameters and ratio Qj/Qr allow any 
value of 9x in the interval (0, it). 

It is noteworthy that, near the critical line Qj — QJ 1 ^, the necessary condition is 
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satisfied if and only if Q x ~ t- Letting Qj/QJ 1 ^ = 1 — 5 we get 
7T > > 1" — 



/2(l + a)[(l + a) 2 + / 1 ]/;\ v 
Ji(J 3 + l + o) J 



3. Recession of the Asymmetric Body (I ± ^ J 2 ) 

In this section we give a precise description of recession, or reversal of precession, of the 
rock'n'roller. The definition is based on observational evidence: for initial conditions 
close to pure rocking motion and such that the local maxima (turning points) of the 
tilting angle are in the range (~ 3/47T, tt), the rates of change 4>(tj) and ipifj) at 
times tj, j = l,...,oo where the angle 0(tj) is a (local) maximum 9x(tj), depend 
on the time tj, contrary to the case of the symmetric body. The functions <j)(tj) 
and 4>{tj) have a quasi-periodic behaviour, undergoing changes of sign that translates 
observationally to alternating reversals of the visible precession angles A<fr(tj) (mod 2ir) 
and Ai[> (tj) (mod27r), where 

A(j)(tj)= I <t>(t)dt, Aip(tj)= / ip(t)dt. 

Jtj — 2 Jtj — 2 

Note that the integration is from t,_ 2 to tj, which accounts for a full period of motion. 
In the dynamical region of interest, Ox(tj) € (~ 3/47r,7r), the critical quantities 
defined in §2.5 allow us to understand the behaviour qualitatively, and to predict 
the occurrence of reversals. 

The key observation from numerical simulations is that, in the asymmetric case, 
the Jellett and Routh quantities, (16) and (17), cease to be conserved, but oscillate 
about mean values. We thus define the Jellett and Routh quantities Qj and Qr to 
be 

Q.,(t) = hsJ 2 + I 3 fuj' 3 , Q R (t) = ^ . (28) 

We have observed that these quantities oscillate about time-averaged values with a 
period that is generally longer than the period of the rocking motion, and that depends 
on the amplitude of the motion. We remark that the motion in the (Qr, Qj)-planc 
is bounded. We will perform a numerical study of this behaviour in connection with 
reversals at the end of next section. The analytical study of this will be the subject 
of forthcoming work. 

The analysis in §2.5 regarding the asymptotic behaviour of <j> near the turning 
points, remains valid if we consider Qj(t) and Qn(t) to be functions of time. In 
particular, as long as the point (Qn(t),Qj(t)) remains within one of the Regions I to 
IV, we can safely conclude that there is no reversal of the system, because the sign 
of <fr at the turning points cannot possibly change. Reversal is due to crossing of the 
system from one region to an adjacent one. In order to observe reversal, we need 
to initialize the system sufficiently close to the boundary of a region in such a way 
that, during the evolution of the motion, the system crosses the boundary. We call 
this a critical crossing. Due to the oscillating nature of Qj(t) and Qr(€) evidenced in 
numerical simulations, if this critical crossing happens then the system will eventually 
cross back to the original region and will continue crossing periodically back and forth 
between the two regions, in a bounded motion within the space (Qr, Qj)- 
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Corresponding to the critical crossings of the two types of critical quantities - 
Q c Jq defined at the turning point near 9 = and Q c j l * defined at that near 8 = w — 
there are two types of reversal. On the one hand, the full angle, 4>{t), has reversals that 
arc related to the critical crossings of Q c Jq- This is due to the fact that, for motion 



close to pure rocking, the main change of (j>(t) from t 



tj-i to t 



tj is typically 



a jump of magnitude about tt when 9 passes the turning point On- The sign of this 
jump depends on which critical region the system is in, and will therefore change when 
reversal occurs. Critical crossings from Region I to Region II or from Region III to 
Region IV correspond to this type of reversal. 

On the other hand, the visible precession angle, A(f>(tj) (mod27r), is due mainly 
to the change of (f>(t) near the turning point Ox- The sign of this change depends 
exclusively on the criticality Q c ]^- Critical crossings from Region IV to Region I or 
from Region II to Region III determine this type of reversal. This reversal corresponds 
to the recession evident in real experiments. A numerical study of the two types of 
reversal will be presented in §6 below. 



4. Asymmetric Body (I ± ^ 7 2 ): The Dynamical Equations 

We now derive the equations for the asymmetric case I\ ^ I 2 - Since the intermediate 
frame (i',j',k') is no longer a principal frame, it is convenient to use a body frame 
k) aligned in the direction of the principal axes. The angular velocity and angular 
momentum are then 

u> = u>ii + uj 2 j + u> 3 k L = hujii + I 2 cu 2 j + h^ 3 k . 

The momentum equations of motion in the body frame are 

vi + lu 2 v 3 - lu 3 v 2 = F± 

v 2 + lu 3 vi - luiv 3 = F 2 (29) 
v 3 + uoiv 2 - uo 2 vi = F 3 
and the angular momentum equations are 
hd>i + (I 3 - I 2 )uj 2 lu 3 = Gi 

I 2 6j 2 + (Ji - I 3 )w 3 wi = G 2 (30) 

I 3 lo 3 + (I 2 - /i)wiw 2 = G 3 

We proceed as in §2, using the constraint (9) to express the velocity v in terms of 
the rotation uj and eliminating the moment G by means of the momentum equations. 
The result may be written 



£0 



KCj 



(31) 



where 




the matrices £ and K are 

X sa 
X -a s X 



—a 




K = 



h 



P 

— i 

-fsa 



-s 2 a\ 



s 2 X 2 



-s 2 ax 
-f 
-.fsx 



s 2 a 2 



-fsa 

-fsx 



(32) 



Precession and recession of the rock'n'roller 



14 



and the vector is 

/ -(g + u\ + w|)asx + {h ~h~ a/)w 2 w 3 \ 
= (g + u\ + Lu%)asa + (I 3 - h+ af)uj\UJ3 , 

\ (h - h)w\u)2 + as(-%wi + 0^2)^3 / 
with x = cosi/j and er = sini/>. Note that neither if nor P w depends explicitly on <j). 
Thus 4> is an ignorable coordinate in the system (31). 



4..1. Special solutions 

We consider pure rocking motion with <f> and tp constant. Then u) = (x#, — a9, 0). 
System (31) implies 

(Ji - h)o- X 6 2 = 

so that a nontrivial solution requires ax = 0. That is, the rocking motion must be 
about one of the principal axes. For a = the system reduces to (12), pure rocking 
about the z-axis. For x = we get the corresponding equation with 7 2 replacing 
I\ and pure rocking about the j-axis. From (13), the ratio of the small amplitude 
oscillations about these principal axes is 

v\ = j h + 1 
V2 V h + 1 ' 

In general, there are no periodic solutions corresponding to the pure rolling motion 
found in the symmetric case. However, if 9 remains zero, we may have spinning about 
the fc-axis, with k vertical. Then lj = (0, 0, ip). The equations (31) reduce to /3W3 = 0, 
confirming that the spin rate is an arbitrary constant. 



4-2. Nonholonomic constraints 



The rock'n'roller is subject to three constraints, one holonomic and two nonholonomic. 
The body must remain in contact with the underlying surface, and the point of contact 
must be momentarily stationary to ensure rolling contact. We can embrace the three 
constraints in the single equation (9), i.e, v = u> X r. We will now express this in 
terms of the space frame. The velocities in the body and space frames, v and V 
respectively, are related by v — Y T V or, explicitly, 

T 

/ Vl \ 

(33) 

where = cos <j> and = sin <p. The matrix T is the product of three rotations, and 
is derived in many standard texts in mechanics; see, for example, [8, 11, 12, 18, 20]. 
We write u> X r = Tu> and u> = where 





c<f>x - cs^a 


C(fiO~ - 


cs^X 


ss<i> 


(•) ■ 


S0X + CC0CT 


-s^a + 


cc^x 






sa 


sx 




c 




r = 





-/ 



/ sx 
sa 
-sa 



(34) 



Now the constraints can be expressed in the form V = TTS 9, relating the velocity 
in the space frame to the time derivatives of the Euler angles. More explicitly, 




hs<j, 
-hc,p 
as 



-ascp -sc^ 
-ass<p -SS0 




(35) 
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It is clear from this form that the constraint on Z is holonomic and may be integrated 
immediately to give Z = 1 — ac = h, the height of the mass centre in terms of the 
tilting angle 9. The constraints on X and Y are nonholonomic. 



5. Lagrangian formulation 

Systems with holonomic constraints can be solved by elimination of redundant 
coordinates or by adding to the Lagrangian a sum of the constraints weighted by 
Lagrange multipliers. When the constraints are nonholonomic, this procedure does 
not apply [11, 20]. We must resist the temptation to substitute (9) into the Lagrangian 
and obtain a Lagrangian that involves only the Euler angles and their derivatives. 
Rather, we must embed the problem in a configuration space of dimension N + M, 
where N is the number of degrees of freedom and M the number of nonholonomic 
constraints. There has been considerable misunderstanding regarding nonholonomic 
constraints; sec [6] for a review. When the constraints are of the form 

g k (q,q,t) = A k (q,t)q = 0, 

that is, where they are linear in the velocities, we can write the equations of motion 
in the form 

d dC dC x-^ dg k 

k 

where /j, k are Lagrange multipliers that can be determined using the constraints. 
In the present case, the configuration space has five dimensions, with coordinates 
(6,<fi,ip,X,Y), the holonomic constraint having been used to eliminate Z. We may 
write the Lagrangian in terms of these coordinates and their time derivatives: 

£ = \{{hX 2 + ho 2 + a 2 s 2 )9 2 + 2{h - h)s X oH 

+ [(ho- 2 + hx 2 )s 2 + I 3 cV + 2hcU + h V> 2 (37) 
+ (X 2 +Y 2 )} -ga(l-c). 

Note that C does not depend on <p. From (35), the nonholonomic constraints are 
g 1 =X-( hs^e - asc^ <j> - sc^ ip) = 

g 2 =Y - (-hc^O - ass^cj)- ss^tp) = (38) 

Although <j) occurs in these expressions, it is absent from the combinations 
Sfc P'kdgk/dq that occur in the equations. This symmetry should imply the existence 
of an invariant quantity in addition to the total energy. 

The Euler-Lagrange equations (36) for X and Y immediately yield 

Mi = -X [i2 = -Y . 

Using the constraints, we may now eliminate the multipliers fj,\ and fj,2 from the 
remaining equations and obtain a system of three equations for 9, (j> and ip. They may 
be written 

M8 + P e (8,0) = (39) 

where = (9, <j>, ip) T and M is a symmetric matrix. The explicit expansion of (39) is 
given in the Appendix (§A.3). Using Mathematica, the system has been shown to be 
completely equivalent to the system (31). 
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In general, we can write the Lagrangian in the form 

£ = £ o (0,0) + e£i(0,V,0) 

where Co is the Lagrangian for the integrable symmetric system and e = (I 2 — h)/Ii 
is the asymmetry parameter. This provides a basis for a perturbation analysis when e 
is small, which will not be undertaken here but will be the subject of future work. 

6. Numerical Experiments 

The numerical integration of the equations is delicate, as there is a singularity of the 
coordinate system when 9 = 0, and significant errors may result from this. To be sure 
of reliable numerical results, we used a routine of eighth-order accuracy, ODE87, coded 
by Vasiliy Govorukhin (http://www.mathworks.com/matlabcentral/), which is a 
realization of the formulae of Prince and Dorman [14]. With this method, invariants 
of the motion remained constant to high accuracy. To further confirm the robustness 
of the numerics, we coded both sets of equations, the system (31) in terms of (9, u>) 
and the system (39) in terms of (9,0), and compared the results. Furthermore, we 
verified the MATLAB coding by an independent coding in Mathematica. Finally, the 
results presented below were checked for convergence by varying the error tolerance. 
We can therefore be confident in the reliability of the numerical results. 

Unless otherwise stated, the numerical values of the parameters are set as follows: 
gravity g = 9.87, unit mass, unit radius, centre of mass off-centering a = 0.05, 
moments of inertia Ii = 0.35 and I 3 = 0.4. Some initial conditions will not be 
varied in the various simulations; these are 9q = 0.957T, 9q = and <po = 0. 

6.1. The consequence of asymmetry 

We first compare the numerical solution of equations (31) for the symmetric case 
h = h and for a case of slight asymmetry. The solutions arc for 1500 time units and 
the initial conditions are, in each case, 9 n = 0.957T, </> = ipo — 0, u>ifl = 0, Ld 2 ,o = 0.001, 
= —0.001. Figure 5 (top left panel) shows the trajectory of the point of contact 
for the symmetric case I\ = I 2 . The azimuthal angle increases regularly and steadily 
for each cycle of rocking motion. This is confirmed by figure 5, top right panel, which 
shows (f> at the points where 9 reaches a maximum. For the solution shown in the 
bottom panels of figure 5, the only difference is a increase of 0.1% in the inertial 
moment about j, so that e = (I 2 — h)/h — 10~ 3 . The bottom left panel shows the 
trajectory of the point of contact: the precession is no longer uniform. The azimuthal 
angle alternately increases and decreases (figure 5 bottom right panel). We see that 
there is recession, with a period much longer than that of the rocking motion. Thus, a 
minute change in the mass distribution of the body, that changes the inertial structure 
slightly and breaks the symmetry 1\ = I 2 , has a dramatic effect on the character of 
the motion. 

6. 2. Stability of rocking motion 

We initiate the motion from a stationary state with u>(0) = and 9(0) = 0.957T. 
Clearly, a symmetric body started in this configuration would execute pure rocking 
motion, passing repeatedly through the equilibrium position, with <f> and ip remaining 
constant (apart from jumps of n due to the coordinate singularity at 9 = 0). For the 
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asymmetric body, the solution depends on the initial angle tp(0) = tpo- As before, we 
assume the asymmetry is slight, with I 2 = l.OOlii. 

The trajectory of the point of contact of the rock'n'roller is shown in figure 6 
for Vo in the set {tt/100, tt/8, tt/4, 3tt/8, 3.9tt/8, tt/2}. All integrations arc for 1000 
time units. We see that the motion precesses through an angle $ = max — 4> m i n that 
depends sensitively on the initial phase V'o- It appears that the relationship 

$ = 7r - 2^o for Vo £ (0: t) 

is satisfied, at least approximately. 

The cases V'o = and Vo = tt/2 correspond to pure rocking about the principal 
axes with moments of inertia I 2 and I\ respectively. Motion close to pure rocking 
about the I\ axis is stable (figure 6(E)) while that starting close to the I 2 axis changes 
dramatically, precessing through almost 180° (figure 6(A)). We recall the classical 
result for free motions of a rigid body with 7i < I 2 < I3, where rotation about the I 2 
axis is unstable whereas rotations about the 1\ and axes are stable. 

In general we expect the trajectory to be dense in the domain of angle <& spanned 
by the solution. However, KAM theory [1, 13] suggests that for exceptional initial 
conditions the solution is periodic. The character of the solution for ?Ao = 37r/8 
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Figure 6. Trajectory of the point of contact in the XY-planc for initial phase 
angle tp £ {tt/100, tt/8, tt/4, 3tt/8, 3.9tt/8, tt/2}. All integrations arc for 1000 
time units. The small circles indicate the starting position in each case. 



appears to be close to a periodic solution (figure 6(D)). Searching in the neighborhood 
of this solution, we found that when tjjQ — 2.9657r/8 the trajectory becomes periodic, 
repeatedly tracing out the same track, some fifteen times in 1000 seconds. Solutions 
of this nature, whose trajectories span a set of measure zero, are a signature of 
integrability. 

6.3. Recession and Criticality 

We now present a numerical study of reversals based on the theory of criticality 
described in §2.5. In all cases, the numerical experiments consist of releasing the 
rock'n'roller at an angle 8q = 0.957T, with 9q — and with an angle ipo = 7r/4 half-way 
between the body's principal axes. System (31) is integrated numerically for 200 time 
units, using an adaptive Mathematica code (stiffness-switching method) with 11th 
order accuracy. The results are insensitive to resolution improvements. We monitor 
energy conservation point-wise and confirm that the relative error is less than 10~ 10 . 
Routh's and Jellctt's quantities are computed in post-processing. 
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Figure 7. (Colour online). Top frame: orbit in the (Qr, Q j)-plane for 200 
time units. Middle frame: Azimuth angle <f>(t). Bottom frame: Visible angle 
4>(t) (mod7r) sampled when 8 > O.50Q. Initial conditions Qr = 0,Qj = 0. For 
full details, see text. 
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Figure 8. (Colour online). Top frame: orbit in the (Qr, Q /)-plane for 200 
time units. Middle frame: Azimuth angle <f>(t). Bottom frame: Visible angle 
<t>(t) (mod n) sampled when 9 > O.50Q- Initial conditions Qjj = 0.0030, Qj = 



^crit w -0.0017. For full details, see text. 
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Figure 9. (Colour online). Top frame: orbit in the (Qr, Q j)-plane for 200 
time units. Middle frame: Azimuth angle <f>(t). Bottom frame: Visible angle 
<t>(t) (mod n) sampled when 8 > O.50o- Initial conditions Qjj = 0.0030, Qj = 
QJ$ « 0.0016. For full details, see text. 
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In figures 7, 8 and 9 all parameters are identical except for the initial conditions. 
In all three cases, 8 = 0.95-7T, 6 = 0, <fr = 0, -00 = 7r /4- Parameter values are 
a = 0.05, g = 9.87, I 3 = 2/5, h = (1 - 5a/2)7 3 , e = (I 2 - h)/h = 10~ 3 . The top 
frame in each case shows the orbit in the (Qr, Qj)-plane for 200 time units (zigzagging 
bounded curve (black)). In all three cases, the orbit starts at a point on one of the 
critical lines and begins moving to the bottom right, subsequently alternating between 
adjacent critical regions. The straight lines from top left to bottom right (red) denote 
the critical line Qj — QJ^, useful for visible angle reversal. The straight lines from 
bottom left to top right (blue) denote the critical line Qj = Q c Jq, useful for full angle 
reversal. Dots correspond to instances when 9(t) is near the turning points: red dots 
denote 9 > 0.99#o (near turning point 9x) and blue dots denote 9 < 0.005#o (near 
turning point 9n)- The middle frame in each case shows the azimuth angle <p{t). Solid 
lines (red) denote instances when 9 > 0.56*o and dashed lines (blue) denote instances 
when 9 < 0.5#o- The bottom frames show the visible (half) angle (f){t) (mod it) sampled 
when 9 > O.50o- 

In figure 7, the initial velocities are </>o = tpo = 0. The motion remains close to 
the centre Qr — Qj — 0, i.e., close to pure rocking (black zigzagging orbit in top 
frame). The system alternates between Regions I and III, spending approximately 
five periods of rocking motion in each region. Consequently, both full-angle reversal 
(middle frame) and visible angle reversal (bottom frame) can be observed. This case 
corresponds exactly to case (C) in figure 6. 

In figure 8, the initial conditions arc right on the critical line Qj — Q c j^- the 
initial velocities are (j>o — —0.002 and ipo = 0.002. The orbit is very similar in shape 
and size to the orbit in the previous case. We observe two critical crossings between 
Region I and Region IV, corresponding to visible angle reversals (bottom frame). 
There is no full angle reversal (middle frame). 

In figure 9, the initial conditions are right on the critical line Qj = Q c ]q'- the 
initial velocities are <po — 0.379 and tpo — 0.378. The orbit differs in shape from the 
ones seen above and its horizontal dimension is three times smaller. We observe eleven 
critical crossings between Region I and Region II, corresponding to full angle reversals 
(middle frame). There is no visible angle reversal (bottom frame). 

We notice in each of the three cases that, regardless of the apparent complexity 
of the orbits, when the system is near one of the turning points (9 = 9x(tj), red dots; 
9 = 0N(tj), blue dots), the points (Qr,Qj) are distributed along a straight line (top 
frame in each case). In figure 10 a plot is shown combining the orbits of the three 
initial conditions used in figures 7, 8 and 9, in order to compare their distribution and 
extent in the plane (Qr,Qj)- 

It is evident from the above that the criticality criterion is a useful description 
of both full angle and visible angle reversals. One just needs to initialize the system 
near a critical line and the dynamics will do the rest. However, we do not yet have an 
explanation for the extent of the orbit, so our method is only descriptive and cannot 
predict, for example, the number of rocking cycles executed in each critical region. 
Forthcoming work should be dedicated to this subject. 

Regarding visible angle reversals, we have found that these cease to be observed 
if the initial (Qr, Qj) is chosen far enough from the origin (keeping all other initial 
conditions fixed). This can be understood from the fact that the asymptotic Laurent 
expansion in (25) is valid only near 9 = n but the maximum attainable 9(t) is bounded, 
from energy conservation, by 8q, which is strictly less than tt. From the analysis given 
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Figure 10. (Colour online). Combined plot of orbits for the three initial 
conditions described in figures 7, 8 and 9. 
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Figure 11. (Colour online). Left frame: high-Qij reversal orbit, predicted by 
the asymptotic theory. Color code as in figure 7. Right frame: visible angle 
4>vaod2iT, sampled at times tj when 9(tj) attains its maximum value 9x{tj). 



in §2.6, and equation (27), it follows that the necessary condition for validity of the 
Laurent expansion becomes Qr < 0.1 for the present choice of parameters and initial 
conditions, where we have used the observational estimate (from figure 10) of 0.001 for 
the orbit extension along the Q^-axis. We have checked that there is indeed reversal 
for Qr = 0.1 (see figure 11, left frame). It is important to mention that at this 
relatively high value of Qr (and correspondingly high angular velocity) the lowest- 
order Laurent asymptotic expansion given in (25) needs to be improved. As a result, 
the simple interpretation of reversals in terms of critical crossings and changes of sign 
of tf> will change slightly. In practice, to observe recession in this limiting case, it is 
necessary to offset slightly the initial condition in the plane (Qr, Qj), to a point above 
the critical line Qj = Q c j^- The resulting orbit remains in Region I so that there is 
no change in sign of <fi(tj). However, the visible precession angle Acf>(tj) (mod27r), 
being determined by an integral in time, can and does have reversals (figure 11, right 
frame) . 
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7. Conclusion 

Experiments show that the recession, or reversal of precession, is a robust feature of 
the motion of the physical rock'n'roller. Analysis has confirmed that for a perfectly 
symmetric body with I\ = I2 this behaviour is impossible. However, even the slightest 
breaking of this incrtial symmetry is sufficient to change the character of the solution 
profoundly, allowing entirely new types of motion. Physical experiments and numerical 
results show that the reversal angle $ = <fi max — <p m in depends sensitively on the 
initial conditions. For motion that is initially close to pure rocking, the angle <f> 
can be controlled by the choice of the initial phase angle ip . A rigourous analytical 
demonstration of this result is outstanding. 

The symmetric equations are integrable, with three invariants: the total energy E, 
Jellett's quantity Qj and Routh's quantity Qr. In the asymmetric case, only one of the 
above three quantities is conserved, namely the total energy. We present an analysis 
of recession based on the existence of critical lines dividing the (Qr, Qj)-plane into 4 
dynamically disjoint regions. We prove that recession is directly related to the lack 
of conservation of Jellett's and Routh's quantities, by identifying individual reversals 
as crossings of the orbit (QR(t),Qj(t)) through the critical lines. The criticality 
criterion allows one to produce a family of initial conditions so that the system will 
exhibit recession. 

In the asymmetric case, there remains an underlying geometric symmetry - 
invariance under change of the azimuthal angle 4> — so it is arguable that another 
dynamical invariant exists. In the realistic case where gravity is present, this additional 
integral (if it exists) remains to be found. Notice that Borisov and Mamaev [2] indicate 
in their Table 1 that the quantity M 2 — 2Kr 2 is conserved in the asymmetric case 
(where M is the angular momentum about the contact point); however this is only 
true in the absence of gravity. 

There is apparently a slow period of the orbit (Qii(t),Qj(t)). This suggests 
multi-scale analysis as an appropriate technique for analysis of this problem. For 
small asymmetry e = (I2 — I\ ) / 1\ the problem may be formulated as a perturbed 
integrable Lagrangian system, and is amenable to standard asymptotic analysis. This 
will be the subject of future work. KAM theory [1, 13] would indicate that certain 
aspects of integrability should apply to the weakly asymmetric rock'n'roller. However, 
the question of the general integrability of the system remains open. 
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Appendix A 

A.l: The basic parameters of the rock'n'roller 

Let us assume that the body consists of homogeneous material of uniform density, and 
that its mass and radius are both unity. We denote by O the co-latitude of the polar 
cap that is removed to construct the rock'n'roller. All the dynamical parameters are 
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determined once this angle is fixed. We define the distance from the geometric centre 
to the centre of the planar face of the body: 

d = cos 

The volume of the body is then 

v = .(l + d-ld* 

The off-set of the mass centre from the geometric centre is 

7T (d 2 {2-d 2 ) -r 



4 V V 

The moments of inertia about the geometric centre are 

/ A + rf _| rf 3 + l d 5 x / _| + l rf+ l rf 3_|, d 5 

73 ~2\ V J' 7l_7r l K 

By means of the parallel axis theorem [18], the moments of inertia about the centre 
of mass are 

h =Iz h=l[- a 2 

For the actual rock'n'roller shown in figure 1, the polar angle is w 53°. Thus d = 0.6, 
giving the (nondimensional) parameter values 

a = 0.085 h = 0.362 I 3 = 0.42 

For our numerical experiments we used the values a = 0.05, 1\ = 0.35 and I 3 = 0.4. 

A. 2: The equations for the symmetric rock'n'roller 

The equations for the symmetric case I\ = I2 were given in §2.1. The details are given 
here. Taking the cross-product of r with the momentum equation (5) gives 

rxi) + rx(flxv) = r X F = r X W - G (A.l) 

Noting that 6 = u>i, the acceleration in i'j'fe'-components is 

* = (/^2 - st ^3> -fa'n sCo 'i) + + CUJ ' 3 ), SUj'i,C(j'i) . 



It follows that 



and 



Moreover, 



rxv = [(s 2 + / 2 )wJ,/ 2 w£ - fsd>' 3 , -fs6j' 2 + s 2 u' 3 ) 

+ (asuj' 2 , —fuj'^su)^ + clu 3 ) , suj' 1 (suj' 2 + cu)' 3 )] 

rx(flxv)r = (r ■ v)fl - (r-fl)v 

= - 2 + cf)(u)' 2 /s)(fu) 2 - suj' 3 , -fu'^sw'i) 



r X W = —gasi' 
Using these expressions in (A.l) we get 
Gi = G° - ( s 2 + f 2 )^ 
G 2 =G° 2 - (f 2 u; 2 - fsu;' 3 ) 
G 3 = G° 3 - (s 2 ^ - f S Lo> 2 ) 
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where, defining the height of the centre of mass as h = 1 — ac, 
G\~— [asto'i — h(fu' 2 — suj' 3 )lu' 2 /s] — gas 

G° 2 = - [-M(SLJ 2 + CUj' 3 ) + hfw[w' 2 /8] 

G 3 = — [sw^(sw 2 + cuu 3 ) — hui'^^] 
We can now substitute for G in (8) to obtain 

[h + (s 2 + fM = -(I 3 n 2 w' 3 - hfl 3 w' 2 ) 



(hn 3 iu[ - h^ 3 ) + g° 2 = 



[h + f}oj' 2 + [-fs}w> 3 



[-fs)uj 2 + [h + s 2 W 3 = 
The first equation immediately gives the evolution of ui[: 

Pi 



Pi 
P 2 
Ps 



(A.2) 



Si. 



1 h + s 2 + f2 

The second and third equations can be written 



h + f 2 

-Is 



-fs 
h + s 2 



2 - 



The matrix is nonsingular, with determinant A = {hh + hs 2 + hf 2 ) and inverse 



m 


1 




~ A _ 



h + s 2 



fs 



fs h+f 2 _ 

The complete system of equations for the angular variables is now obtained: 
ui't , <p = u> 2 /s, ip = u> 3 — (c/s)u>' 2 . 



■ 1 



J l ! 

Si 



UJ 9 



S2, 



This system provides six equations for the six variables {9, (f), tjj, to[ , u)' 2 , ui' 3 }. 
A.3: The Euler- Lagrange equations 

The Lagrange equations arising from (38) may be written 

M0 + Pg(O,b) = 

where = {9, <j>, ip) T . The symmetric matrix M is defined as 
" hx 2 + ho 2 +a 2 s 2 + h 2 (h-I 2 )so-x 



(A.3) 



M = 



{h - h)sax 




{ho 2 + hx 2 + a 2 )s 2 + he 2 
he + as 2 





he + as 2 
h + s 2 



and the vector Pq = (Pg, P^, P^) has components: 

P e = [as]6 2 + [-{ho 2 + I 2X 2 - h)sc + has}<j> 2 + [-{h - h)2o X ]H 

+ [{h - h)s{x 2 - o 2 ) + I 3 s + hs]<j>i> + gas , 
P = [{h - h)co X ]9 2 + [{ho 2 + I 2X 2 - h)2sc + {2ac - l)as]64> 

+ [{h - h)s{ X 2 - o 2 ) - hs + asc]6i> + [{h - h)2s 2 a X ]H , 
P^ = [{h - h)o X ]0 2 + [-{h - h)s 2 o X ]4> 2 

+ [-{{h - h){x 2 - o 2 ) + h)s + {2ac - l)s\H + [sc]6^ . 

Now (A.3) may be solved for {0{t),6{t)). It has been confirmed, using Mathematica, 
that the system (39) is completely equivalent to the system (31). 
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